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Abstract 

We study the one-dimensional Dirac equation in the framework of a position dependent mass under the 
action of a Woods-Saxon external potential. We find that constraining appropriately the mass function it 
is possible to obtain a solution of the problem in terms of the hypergeometric function. The mass function 
for which this turns out to be possible is continuous. In particular we study the scattering problem and 
derive exact expressions for the reflection and transmission coefficients which are compared to those of the 
constant mass case. For the very same mass function the bound state problem is also solved, providing a 
transcendental equation for the energy eigenvalues which is solved numerically. 



PACS numbers: 03.65.-w; 03.65.Ge; 12.39.Fd 



I. INTRODUCTION 



In recent years, the study of several quantum mechanical systems within the framework of an 
effective position-dependent mass (PDM) has received increasing attention in the literature. Posi- 
tion dependent effective masses enter, for example, in the dynamics of electrons in semiconductor 
hetero-structures [I], and when describing the properties of hetero-j unctions and quantum dots [2]. 
In non-relativistic quantum mechanics when the mass becomes dependent on the position coordi- 
nate the mass and momentum operators no longer commute, thereby making the generalization of 
the non-relativistic Hamiltonian (kinetic energy operator) to the PDM case highly non trivial [31 2] , 
because of the ambiguities in the choice of a correct ordering of mass and momentum operator [5] . 
Another important issue is that of Galilean invariance [6]. 

The investigation of relativistic effects is of course important in those systems containing heavy 
atoms or heavy ion doping [7j. Therefore for these type of materials the investigations of the 
properties of the Dirac equation in circumstances where the mass becomes a function of the position 
is certainly of great interest. In addition the problems posed by the ambiguities of the mass and 
momentum operator ordering are absent in the Dirac equation. An effort in this direction has 
been reported in some recent literature [7Ml5|. For example the authors of ref. [7] have reported 
an interesting numerical investigation of the scattering problem for the three-dimensional Dirac 
hamiltonian within a position dependent mass with a costant asymptotic limit, studying the energy 
resonance structure. In Ref. |16j the scattering problem is solved for a smooth potential and a mass 
step but in the non-relativistic regime. The authors of ref. |17j reported an approximated solution 
of the one-dimensional Dirac equation with a position dependent mass for the generalized Hulthen 
potential. To the best of our knowledge few attempts have been reported that study the Dirac 
equation in an external potential with position dependent effective masses. In ref. [18] the author 
studies Dirac equation in 3+1 dimension in the Coulomb field and with a spherically symmetric 
singular mass distribution. In Ref [10J the author reports an exact solution for the Dirac equation 
with central potential and mass distribution both inversely proportional to the distance from the 
center. 

It is worth pointing out that graphene (single atomic layer of graphite), a recently discovered 
material [19, 20J which is receiving a lot of attention, exhibits several properties whose explanation 
involve the Dirac equation for massless fermions. For a comprehensive review see for example [21] 
Recent reports studying these effects [221 [23] attest the use of the Dirac equation in explaining the 
properties of single-layer graphene . 
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The present work is an attempt in the same direction of ref. |1CH [T8] . We report an exact solution 
of the one-dimensional Dirac equation in the position-dependent mass formalism for a particle in 
the Woods-Saxon (WS) potential. Our approach is based on that of Ref. [21] where the author 
solves the one-dimensional Dirac equation in the Woods Saxon potential for the ordinary constant 
mass case. Our method consist in requiring, within the effective position dependent mass, that 
the second order equation still be exactly solvable by the hypergeometric function. This is done 
imposing restrictive conditions on the mass function which lead to a first order differential equation 
which provides the explicit mass function. 

We would also like to stress that our new analytical exact solution of the position dependent 
mass Dirac equation in the Woods-Saxon potential will prove certainly of use in further studies 
of effective mass models. Other issues could for example be addressed that go beyond the scope 
of the present work: for example it would certainly be of interest to study in the detail the Klein 
paradox, as well as the issue of zero momentum resonances that support a bound state at E = — m, 
i.e. super-critical states, in the framework of effective masses. We also note that the study of the 
transmission coefficient in the case of two dimensional Dirac equation for massless fermions has 
been used already to describe the electrical properties of graphene and in particular the possibility 
to observe the Klein paradox phenomena |25j in this material. In addition the authors of [23] 
discuss the case of massless electrons that cross a square barrier region where they are instead 
massive, a situation that can simulate a n — p — n junction in a graphene nano-transistor. We 
believe that our exact solution derived for the Woods-Saxon potential with an effective position 
dependent mass, in the limit aL S> 1 the WS potential barrier reduces to a square barrier, may 
prove useful to describe such real physical system. Further investigation in this direction is needed 
but is beyond the scope of the present work. 

The plan of the paper is as follows. In Section II, we summarize the basic equations of the 
problem. In Section pTIj we solve the effective-mass Dirac equation and provide the mass function for 
the problem. We study the scattering problem for the potential barrier and deduce the transmission 
and reflection coefficients by studying the asymptotic behavior of the wave function when x —> ±oo 
and the match at x = 0. In section [TV] we also address the bound states of the problem by turning 
the Woods-Saxon potential barrier into a Woods-Saxon potential well. We discuss several numerical 
examples providing the eigenvalues and wave functions corresponding to particular choices of the 
parameters. Finally we summarize our results and present our conclusions in Section W\ 



3 



II. THE BASIC EQUATIONS 



We recall the free Dirac equation (using natural units, h = c = 1} 

,. d 



m *f?(x) = 0. 



(1) 



and [i = 0,1,2,3 is a space-time index. Considering instead the case of one space dimension 
it is possible to choose the gamma matrices j x and 7 of dimension two and it is customary 
to set them respectively to the Pauli matrices ia x and a z [9]. Considering a charge particle 
minimally coupled to an electromagnetic potential, in the absence of the space component of a 
vector potential, and setting V(x) = eAo(x) the one-dimensional Dirac equation for a stationary 
state ^(x,t) = e~ lEt ij){x) becomes: 

d 



Gx ~Ax ~ ( E - V ( x )) a z + m 



if){x) 



0. 



(2) 



Decomposing the Dirac spinor ip{x) into upper (u±) and lower (112) components: ip = [ ~* J, gives 



the coupled equations: 



Ui(x) = — [m + E — V(x)] U2(x) , 
u' 2 (x) = — [m — E + V{x)\ u\(x) , 



(3a) 
(3b) 



It turns out to be convenient to define two auxiliary components <j>(x) and x( x ) i n terms of u\{x) 
and U2(x) as in [9]: 



4>(x) = u\{x) + lU2(x) 
X(x) = ui(x) - iu 2 {x) 



(4a) 
(4b) 



Using the above definitions and Eqs. ( 3apb ) we find the first order coupled equations for the 
components (p(x) and x(x): 



4>(x) = —imx(x) + i [E — V(x)] <f>(x) 
x'(x) = -\-im(f>(x) — i[E — V(x)] x( x ) 



(5a) 
(5b) 



which give the second order equations: 



4>"{x) + [{E - V{x)) 2 -m 2 + iV'(x)] <t>(x) = -im! x(x) , 
X"(x) + [(E - V(x)) 2 - m 2 - iV'(x)} x( x ) = +im' <p{x) , 



(6a) 
(6b) 
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where we have taken into account the fact that the mass may depend on the position coordinate 



and prime denotes the derivative with respect to x. Eliminating x( x ) using Eq. (5a), we obtain for 
4>{x) the second order equation: 



in 



fix) - —ct>'(x) + 



m 



(E - V(x)) 2 -m 2 + iV'(x) + i—(E- V(x)) 

m 



0. 



(7) 



Solving this second order differential equation one can obtain the x( x ) component via Eq. (5a) 
and then reconstruct the upper u\{x) and lower U2{x) components of the complete spinor solution 
ip{x). Eq. ([7]) reduces to the one studied in [24] if m' = 0, i.e. if the mass reduces to a constant. 
Thus we see that keeping a position dependence in the mass introduces two new terms: one which 
multiplies (j)'(x) while the other enters the <j>[x) term. These new terms must be appropriately 
constrained in order to be able to solve the equation in terms of the Hypergeometric function. 

Let us make a final remark before to discuss the details of the computations. The attentive 
reader might wonder what would happen if one were to derive a second order equation for the 



x(x) eliminating instead the <fi(x) component and then computing it via Eq. 5b The second order 
equation for the x( x ) component turns out to be: 



m 



x"{x) x'{x) + 



m 



m 



(E - V(x)y - m 2 - iV'(x) -i — (E- V{x)) 



rn 



X(x) = 0. 



(8) 



It is easily checked that Eq. [8] can be obtained from Eq. [7] using the map E — > —E and V(x) — > 
—V(x). This can be interpreted as the negative energy solution corresponding to the charge 
conjugate particle (antiparticle). Since V{x) is the temporal component of a four-potential the 
change V{x) — > — V(x) amounts to reversing the charge of the particle. Indeed with E — > —E 
and V(x) — > —V(x) we have x ~ * 4> > an d by using Eq. 5a, (j) — >■ — x and by using the inverse 
of Eqs. fl4a|4b"l ) we have in turn u\ — > iu2 and 112 —> iu\ which amounts, up to inessential phase 
factors, to the charge conjugation symmetry of the Dirac equation 



III. EFFECTIVE-MASS DIRAC SCATTERING PROBLEM 



The form of the Woods-Saxon potential (illustrated in Fig. [I]) is given by (see also ref. |24j ) : 

9{-x) 9{x) 



V(x) = W 



+ 



(9) 



e -a(x+L) _|_ ^ e a(x-L) _|_ -y 

where W is a positive parameter in the scattering problem (potential barrier) and negative in the 
bound state problem (potential well); a and L are two real and positive parameters. 
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A. Solution in the negative region (x < 0). 



Using the variable y = —e~ a ^- x+L ^ and the transformation eft = y^(l — y)~ X f(y) used in (23j we 
will show that it is possible to obtain an exact solution in the form of an hypergeometric function 
by imposing appropriate constraints on the mass function. With the above transformation and 
using Eq. (8), Eq. ([7]) becomes: 



2/(1 - y) 



d 2 f(y) 
dy 2 



+ 



Tfl 

1 + 2 M - y(l + 2 M - 2A) - - y(l - y) 

m 



df(y) 

dy 



M 2 (l-y) 2 + A(l + A)y 2 + 



2/(1 - y) 

+ [{E 2 - m 2 )(l - yf + W 2 - 2EW{1 - y) - iayW] 

W 



a- 
til 



m 



2/(1 - y) 



A + -\E 



f(y) = o. 



(io) 



y 1-2/ ay V 1 - 2/, 

Here and in the following the dot indicates derivation with respect to the transformed variable 
(m = dm/dy). In order to keep the structure of the hypergeometric differential equation we 
impose the following condition on this term: 



m 



m 



2/(1 - y) = a + (3y 



(11) 



which has the following solution (mo integration constant): 

|y_ i\a+P 



m(y) = m 



\y\ 



(12) 



With this choice of mass function Eq. 10 becomes: 



V^-y)^^- + [l + 2^-y(l + 2^-2X) + (a + f3y)] dJi!l) 



dy 2 



+ <^A(l + 2/i) + 



1 



2/(1 - y) 



dy 

/i 2 (l-y) 2 + A(l + A)y 2 + 



+ -j [(£ 2 - m 2 )(l - y) 2 + W 2 - 2EW(l - y) - iayW] 



n(l-y) + \y + -(E(l-y)-W) 



f(y) = o (13) 



In order that Eq. 13 keeps the structure of the hypergeometric differential equation as in the 
m = const case we may impose the following conditions on the mass function: 



lim m(y) = mo 

y— >— oo 

2/i \2 2 2 

™ (i - y) = m oy 



(14) 
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From these conditions we fix completely the two parameters to a = — 1 and f3 = so that the mass 
function (in the y < region ) becomes: 



m(y) = m 



(15) 



The most general condition that we can impose on the term multiplying to l/[y(l — y)] in order that 
the equation be that of the hypergeometric function is that it be equal to a constant 7. Therefore 
we get three equations: 



E 2 W 2 



EW 



fj 2 -\ 5- + 2 — ^ fi — -(E— W) = 



a 2 a 2 



/;- ./:ir ,w-e , 2 



2 + 2 " - 



2^ - (A - (j) = 7 

-7 



(16a) 
(16b) 
(16c) 



From Eq. (16a), it is possible to solve for \i while A is found summing Eq. (16b) and Eq. (16c). We 



obtain finally: 



A = i 



W 2 



m. 







a* 

(E-W) 



a 



7 = v 2 - fj 2 - A(A + 1) 



(17a) 
(17b) 
(17c) 



having defined v = ik/a where k 2 = E 2 — m^. Our Eq. (13) becomes the differential equation of 
the hypergeometric function 

y(l-y) d ^ + [2fi-(l + 2 f ,-2\)y] d ^-( f ,-\-v)( i ,-\ + v)f = 0, (18) 

and the general solution is (with D\ and D2 constants): 



f(y) =D 12 F l (LL-v-\,n + v-\;2 f i;y)+D 2 y 1 - 2 >* 2 F 1 (l-LL-v-\,l- t i + v-\;2-2fi;y). (19) 
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B. Solution in the positive region (x > 0) 



Let us study the other region in which x > 0, from Eq. fusing the variable 1/ z = 1 + e a ^ x L ^ 
and the transformation <p = z~ u (l — z)~ p g(z), we obtain: 

d 2 g{z) \ , . m ( .1 dg(z) 



zl-z) 



dz 2 



+ 



Tfl 

1 - 2uj - z(2 - 2p - 2uj) z(l- z) 

m 



dz 



m 



+ \ 2uj{l-p) z(l-z) 

m 



+ 



+ 
+ 



1 

1 



z 1 — z 

uj{u + !)(!- zf + p(l + p)z 2 - w(l - z ) + zp- (2z 2 )p 



m 



(E - WzY - m(zY - iaWzCl - z) - ia(E - Wz)—z(\ - z) 

m 



g{z) = 
(20) 



Following a line of thought similar to the one outlined in the previous subsection we obtain the 
mass function 



m(z) = mo(l — z) . 



(21) 



The final result for the parameters cj, p and 5 (the coefficient introduced requiring that the term 
multiplying the factor l/[z(l — z)] be a costant) is found to be: 



{E — W) 
p = -i = p 



uj = i 



E 2 



k 

i — = v 
a 



5 = A 2 - p 2 - v 2 - 2v 



(22a) 

(22b) 
(22c) 



The differential equation of the function g{z) in Eq. (20) becomes: 
d 2 g(z 



z(l-z)- 



+ [l-2u- z(l - 2p-2u)\ 



dg(z) 



dz 2 ' L /J dz 

with the general solution (d± and di constants): 



(-p-L>- A) (-p-u + X) g(z) = 0, (23) 



g(z) = d\ 2 Fi(-p - v - A, -p - v + A; 1 - 2v- z) + d 2 z lv 2 F 1 (-p + v - A, -p + v + A; 1 + 2u; z) . (24) 

Let us briefly comment on the change of variables chosen in the negative region x £ [— oo,0]: 
y = — e - a ( x + L ) anc | i n the positive region x 6 [0, +oo]: 1/z = 1 + e a ( x - L ), This implies clearly 
that y £ [-oo,-e" aL ] and z E [(1 + e~ ai ) _1 ,0] (which reduce to y £ [-00, 0] and z £ [1,0] 
with the assumption aL 3> 1). The reader might worry that the choice in the negative region 
is not appropriate because then the Hypergeometric series 2-^1 (a, b, c; y) would not be convergent 



as its radius of convergence is \y\ < 1. However we would like to stress that (in the negative 
region) we use the explicit form of the Hypergeometric series only in vicinity of x — > CP (|y| <C 1) 
(in the scattering problem) remaining well within the radius of convergence (|y| < 1). When 
obtaining the asymptotic expression at x — > —oo (y —> — oo) we use the asymptotic expansion of 



the Hypergeometric series given in Eq 27 In the positive region the expansion in x — > + (z — > 1 ) 
to avoid problems of convergence we use the continuation identity of the Hypergeometric function 



given in Eq. 36 We can say that our exact solution does not rely at all on the Hypergeometric 
series but rather on the properties of the differential equation satisfied by 2-^1(0, b, c; y). If one needs 
to use the Hypergeometric function outside the radius of convergence of the series it is necessary 
to use adequate analytic continuation identities [26] . We expect in any case that using other 
transformations like for example y = e ax and z = e~ ax would in the end not alter our conclusions. 
Similar considerations apply as well to the change of variables chosen in the discussion of the bound 
state problem (see section IV, subsections A and B). 



Mass function 



From conditions given in Eqs. fll5|21| ), we have obtained the mass function in the two regions as 
follows in terms of the x variable: 



m(x) = mo 



„—a(x+L) p a(x—L) 

■ e(-x) + 



e ~a(x+L) _|_ I 



e a(x~L) _|_ I 



Q(x) 



(25) 



Even if the conditions on the parameters a and (5 are different in the two regions we obtain a 
continous function, infact the limit for x — > + and x — > CP is the same; the value in x = is 
m(0) = mo e~ aL I (e~ aL + 1) < mo as given in Fig.[l} This mass function has the desired asymptotic 
behavior at x — > ±00 where it assumes the desired constant value mo- The reader might be worried 
that derivatives of the 0-functions might introduce singularities in the problem, since in Eq. 
there appear terms with first derivatives of the mass function. However it is quite straightforward 
to check that, due to the continuity of the mass function at x = 0, such 5- function contributions 
do cancel exactly. 

An interesting and relevant point to address comes from the fact that the mass function is of 
the same type of the vector potential (the Woods-Saxon potential). Indeed the mass function that 



we have derived (Eq. 25 can be written in the following form : 



m(x) = mo — r yV(x) 
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FIG. 1: Left, plot of the Woods Saxon potential (W = 1). Right, plot of the position dependent mass 
function ttiq = 1. The parameters are: L = 2 and a = 10 (solid line), a = 5 (dashed line) and a = 3 
(dot-dashed line) . In the limit aL > 1 the Woods-Saxon potential reduces to a smooth barrier approaching 
a square barrier. 

where 7 = uiq/W. Therefore the position dependent mass problem could be looked at as that of a 
particle of constant mass mo coupling to both a vector potential (V(x)) and to a scalar potential 
S(x) = —'jV(x). When discussing the barrier (W > 0) it turns out that 7 > and therefore since 
S(x) + V(x) = (1 — 7)V(cc) there might be regions of the parameter space, mo = W (or mo ~ W) 
where the system is endowed with either an exact (or approximate) pseudo-spin symmetry which 
is denned by the condition S(x) + V(x) = constant. Such symmetry, the near equality of an 
attractive scalar potential with a repulsive vector potential is well know in the literature |27| [28] 
of the Dirac equation and has been proved very useful in describing the motion of nucleons in the 
relativistic mean fields resulting form nucleon-meson interactions, nucleon-nucleon Skyrme-type 
interactions and QCD sum rules. 

Further investigation of the possible consequences of this symmetry for the system under con- 
sideration goes however beyond the scope of the present work. 

D. Asymptotic expressions and boundary conditions of the scattering problem 



In the negative region from Eq. ( 19 ) we have 



b L (y) = D l y^{l-y)- x 2 F l {^-u-\^ + v-\ ] 2 l i-y) 

+ D 2 y 1 ^ (1 - y)~ x 2 Fx(l - fx - v - A, 1 - M + v - A; 2 - 2fj,;y) . (26) 
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We can derive the asymptotic expression as x — > — oo (y —> — oo) by using the following formula 
for the asymptotic behaviour of the hypergeometric function [26J 

^ / , x T(c)T(b — a) , . „ T(c)T(a — b) , . h 

^^y) = f^fF^y <-»r + fUfW ( -*> • (27) 

and obtain: 

(^l(x) ~ G e - tki - x+L) + i? e ^+ L ) (28) 

where 

G = D x A e i7rfi -D 2 C eT™* (29a) 

H = D\B e i7Tl * - D 2 D e~ iirfi (29b) 



and A, B, C, D are given by: 



A - r(2 M )r(2„) , &) 

_ r(2/x)r(-2^) 

5 " r( M - «, - A)r( M - „ + A) (30b) 

C = T{2 - 2/i)F(2l/) (30c) 

r(i-/i + i/-A)r(i-/i + i/ + A) { ' 

D - F(2 - 2p)T(-2u) 



Similarly we can derive the asymptotic form of lower component xi x ) from Eq. (5a): 



lim X (x) L = G e ~*H*+L) + H e ik(x+L) (31) 



m m 



Similarly for the solution in the positive region we have from Eq. (24): 



<f> R {z) = d\ z~ v (1 - z)~ p 2 Fi(-p - v — \, —p — v + \; 1 — 2v; z) 

+ d 2 z v (1 - z)- p 2 Fi(-p + v - A, -p + v + A; 1 + 2u; z) . (32) 

Now we recall that z — > when x — > oo and imposing the boundary condition of the scattering 
problem that in the {x > region) we only have a wave travelling to the right (only the transmitted 
wave) we find: 

lim <p R {x) = d x e lk{x - L) (33) 

x— s>+oo 



and Xr( x ) is found again through Eq. (5a) in terms of 4> R : 



lim XR (x) = d 1 i -^±e^-V (34) 

X^ + OO niQ 
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E. Match of solution at x — 



So far we have derived asymptotic expressions at x — > ±00 for the wave function of the scattering 
problem in the negative region, x < 0, (4>l), and in the positive region, x > 0, (4>r) which depend 
respectively on two, D\ and D 2 , and one, d%, unknown constants. In order to have a physical 
solution of the scattering problem one needs to match the two solutions 4>l arid 4>r at x = 0. This 
is done by imposing the continuity of the wave function and of its derivative at x = which gives 
two conditions and two of the three unknown constants can be expressed in terms of the one left 
out as the ordinary normalization constant. 

We need to find the behavior of the function 4>l(%) and (/)r(x) in the vicinity of x = 0. 

As x — > we have \y\ ~ e~ aL <C 1 as our only assumption throughout the paper is that aL S> 1. 



Thus (1 — y) ~ 1 and from Eq. 26 we obtain for 4>L(y): 



<f> L (x) ~ D^-e-^+^Y + D 2 {-e- a{x+L) ) 1 -^ 

having evaluated the two hypergeometric functions to unity as their argument vanishes, the above 
equation can be put as: 

<t> L (x) ~ D x e in ^ e~ a ^ x+L) - D 2 e ~ a{x+L) e a ^ x+L) (35) 

In order to extract the behavior of 4>r(x) as x — > we use the continuation identity of the hyper- 
geometric function |26j 

2 F 1 (a,b;c;z) = (36) 
T(c)T(c- a-b) 



T(c-a)T{c-b) 



Fna, b; a + b — c + 1; 1 — z) 



+ (1 ~ z) c - {a+b) m 5\rL~ C) 2F l (c-a,c-b;c-a-b+l;l-z). 

r(o)r(6) 

Proceeding similarly to the case x < we find (for x > 0), that as x — > + , z — > 1, and 1 — z 



e a(x ( we assum e aL S> 1) and recalling that /x = p, from Eq. 32 with d 2 = we find:, 

0(aO* ~ d i M e- a i i i»- L ) +d 1 N e (A*+i)«C*-&) ( 37 ) 

where: 

m - r(i-2,)r(i + 2^) 

_ r(i-2^)r(-i-2 M ) 
* - r(-/x - 1/ - A)r(-/i — f + a) (38b) 
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The match of the two solutions </>l(x) and <Pr{x) is done by imposing the continuity of the wave 
function and of its derivative at x = which gives: 



Di ^ e~^ aL + (1 - n)a D 2 e~ i7Tfl e ~^^ aL = di \-fia M e^ aL + N (p + l)a e - (/x+1)aL 



and solving: 



1 - 2/i 



[M (1 - 2/x) e 2 ^ aL + 2iVe- ai ] 
[iV(2 M + l)e~ 2lxaL ] 



(39a) 
(39b) 



1 - 2fi 

We would like to remark that when solving the Dirac equation the continuity condition at a given 
boundary (x = in our case) should be imposed by requiring the match of both the upper and lower 
spinor components (u\(x) and 1*2(2;)). In our second order approach based on the introduction of 
the auxiliary components <j)(x) and x( x )i the derivative of one of the two (<j)') is connected to the 



other (x) because of Eq. 5a In turns both the initial upper u\ and lower u 2 components can be 



expressed in terms of cj> and <fi' . Indeed solving Eq. 4a and Eq. 4b and using Eq. 5a one finds for 
u\ and U2- 



ui(x) 
u 2 (x) 



1 + 
1 - 



E - V(x) 
m 

E - V{x) 
m 



<f>(x) + —<t>'{x) 
m 

4>{x) 4>'{x) 



which shows how the matching of the wave function <fi(x) and of its derivative <f>'{x) is totally 
equivalent to requiring the continuity of 1*1(2;) and 1*2(2;). We also remind the reader that this 
method of matching the solution in x = has been used also in ref. [23]. The above consideration 



applies as well to subsection C of section IV 



F. Probability current density, reflection and transmission coefficients 

The reader might wonder whether the Dirac equation with a position dependent mass still has 
a conserved current. It is well know that a continuity equation for the current is related to the 
conservation of probability, or unitarity. It is quite straightforward to show that a coordinate 
dependence of the mass does not bring in any change in the derivation of the conserved current. 
This is related to the fact the mass multiplies the spinor wave function and in deriving the conserved 
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current such terms simply drop out as in the constant mass case. The probability current density 
is given by: 

J(x) = ij;(x)'y x ij>(x) = i [u*(x)u2(x) — U2(x)u±(x)] 
which can be also given in terms of the auxiliary functions <f> e x : 

j(x) = 1 -[\m\ 2 -\x(x)\ 2 ] 

With the asymptotic form of the wave function we can compute the left (x < 0) Jl, and the right 
(x > 0) Jr current density: 

Jh = Jinc — Jrefl = \H\ 2 — ^ ~ ~~ I^P ~~ " 2 " ' (^0) 

m Q m 

Jr = J trans = \dl\ 2 5 , (41) 

and we can define the transmission and reflection coefficients: 

rp _ Jtrans _ 1^1 1 / 

J ~ \H\ 2 1 ' 

R = Jrefl (E + k) \_Gl_ 

J mc (E-k)\H\* {6) 
and from the current conservation d x J{x) = it follows that f_°° d x J(x) = [J(x)]~^ = Jr — Jl = 
0, and therefore Jl = Jr from which we have the unitarity condition R + T = 1. The transmission 
and reflection coefficients can then be written as: 

|l-2/i| 2 



MB(1 - 2fi) e 2 » aL + N [B (2 — 2/j,) e (- 2 **+i)a£ _ p ^ + 1) e -^ aL ] | 2 



(44a) 



E + k\A[M(l-2fi)e 2 '" lL + N(2-2^)e- aL ]-C[N{l + 2fi)e- 2 ' iaL }\ 2 
~ E-k \B[M (1 - 2[i) e 2 ^ aL + N (2 - 2/z) e~ aL ] - D[N(l + 2//) e- 2 >" aL }\ 2 1 " 

Fig. [2] shows the transmission coefficient in the constant mass case (left plots) and for the position 

dependent mass case (right plot) for two choices of the parameters (a, L) and in the so called Klein 

range m < E < W — m. We note that in the PDM case we still observe the transmission resonances 

found for constant mass. We observe that, while for m = itlq when E — > W — m, T — > 0, in the 

PDM case T — > 1. This is an important fact wortwhile to be pointed out. In Fig. [3] we plot the 

transmission coefficient as a function of the barrier height W. We note that as opposed to the case 

of constant mass where T = for E — m < W < E + m, in the position dependent mass case, T(W) 

always oscillates and does not go to zero in the interval E — m < W < E + m. We note also that 

we have verified our numerical calculations of the constant mass case with those of ref. [24J finding 

complete agreement. Finally, we have numerically checked the validity of the unitarity condition 

R + T = 1. 
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FIG. 2: Left, Transmission coefficient in the constant mass case. Right: plot of the transmission coefficient 
in the position dependent mass case. Parameters: a = 5, L = 10. Upper plots are for W = 1.2 and lower 
plots for W = 4.2. 
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FIG. 3: Transmission coefficient as a function of the barrier height W. Left: in the constant mass case. 
Right: the position dependent mass case. Parameters: a = 5, L = 10, m = 0.4 and E = 2m . 

IV. EFFECTIVE-MASS DIRAC EQUATION, BOUND STATES 



Let us study the bound states for the particle with position dependent mass. In order to do 
this we take the W-S potential with W — ^ —W '. 
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A. Negative region 



In the study of the discrete spectrum it is convenient to use a different variable. Now we choose 
y^ 1 = 1 + e ~ a<yX+L \ with d/dx = ay(l — y)d/dy, V(y) = —Wy and m(y) = mo (1 — y) and using 
the parametric transformation <j> = y a (l — y) e h(y), we obtain from Eq. ^ 



2/(1 - y) 
l 

+ 



d 2 h(y) 
dy 2 



+ [1 + 2a - y(l + 2a + 2e)] 



dh(y) 
dy 2 



+ 



-[(E + Wy) 2 - mg(l - y) 2 - m^y(l - y) - iay{E + Wy)] h(y) 



a 2 y(l - y) 

+ * J g(o- - 1)(1 - y) 2 + e(e - l)y 2 + a(l - y) 2 \h(y) + {-2ea - e)h(y) = (45) 

y(i - y) 

The most general condition that we can impose on the term multiplying l/[y(l — y)\ in order that 
the equation be that of the hypergeometric function is that it be equal to a constant Q. Therefore 
we get three equations: 



a 2 + 



E 2 ^ 

,2 



+ 2- 



EW 



a 

.W + E 



a 



a 2 + e 2 -e + 



W 2 



2a l 





C 

-c 



(46a) 
(46b) 
(46c) 



From Eq. (46a) we can solve for a, while summing Eq. (46b) and Eq. (46c) we obtain the equation 
for e 



a 



Vm 2 ~ E 2 



a 



(E + Wf 



.W + E 2 

i e + e 2 = 



(47a) 
(47b) 



a* a 

So a = v. A solution of the second one is e = —i(E + W)/a, which is the same of \i in the 
scattering problem but with the replacement W — > —W. In determining the parameters a and b 



(c = 1 + 2v) of the hypergeomtric equation we make use of the relation in Eq. (46c), where we 
define A = i \J (W 2 — mfy/a 2 and finally obtain: 

VO- ~ y)^f- + [1 + 2^ - (1 + 26 + 2^)y]^^ - (6 + v + A)(e + v - A)/(y) = 0. (48) 



dy 2 



dy 2 
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B. Positive region 



We consider the same substitution of scattering states for the variable x in the positive region. 
So we use the variable 1/2 = 1 + e a ( x ~ L \ and the transformation (f> = z~ T (l — z)~ n g(z) 

^d 2 q(z) r , .^dq(z) 

z(l - z)— fY- + [1 - 2r - z{l - 2 V - 2r)]- b 



dz 2 dz 

1 

a 2 z(l — z) 



+ ^5— r[(E + Wzf - mg(l - z) 2 + iaWzCl - z) + iaz(£ + Wz)]g(z) 

n Z ~ I 1 



+ 7- - . [t(t + 1)(1 - z) 2 + 77(77 + l)z 2 - t(1 - z) + r]z - 2r]z 2 - tz{1 - z) + 77Z 2 ]5(z) 
z(l zj 

+ (2r - 2T7])g(z) = (49) 
Following a line of thought similar to that outlined in the previous subsection (x < 0) we obtain 



t 2 = —(E 2 — mfyja 2 = o —v and 77 = —i(E + W)/a = e so that Eq. 49 reduces to: 



z(l - ^^fi + [1 - 2u - (1 - 2e - 2i/) -(-e-i/ + A)(-e - 1/ - A)/(*) = 0. (50) 



C. Bound state wave function and match at x — 

We note that the wave function in the x > region can be obtained from that of the x < 
region simply letting v — > —v and e — )• — e. The general solutions to Eqs. ( 48|50 ) are: 



h{y) = A' 2F 1 {e + v + \,e + v-\;l + 2v;y) + B'y- 2u 2 F 1 (e-v + \,e-iJ- X;l-2u;y), 
g (z) = C 2 Fi(-e - v + A, -e - v - A; 1 - 2u; z) + D' z 2u 2 F\ (-e - 1/ + A, — e — i/ — A; 1 + 2z/; z) . 

Recall the parametric transformation for 0L 5 i?: </>_r = z _I/ (l — z)~ e g(z) and <^>£ = — y) e h(y) 
and that in the limit of x — > ±00 the variable y — > as well as z — > 0. Therefore imposing the 
boundary condition of a bound state (vanishing wave function at infinity) we obtain B' = C = 
and we are left with: 

Mv) = A' y^l-yY 2 F l (e + v + \,e + v-\,l + 2v;y) 
<P R (z) = D' z v {\ - z)~ e 2 Fx(-e + v + A, -e + v - A, 1 + 2u; z) 



With the help of the continuation formula of the Hypergeometric function [26] we can extract the 
behavior of the solution in the vicinity of x = (recall that for x — > 0, y, z — > 1 and 1 — y ~ e- a i x ~ L ) 
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while 1 - z » e< x ~ L ">) 



<Pl{x) 
A' 

4>r{x, 

D' 



r(l + 2v)Y(l - 2e) 



r(l - e + 1/ - A)r(l - e + z/ + A) 
r(l + 2z^)r(l + 2e) 



-ea(xH 



+L) + r(l + 2^)r(-l + 2e) p _ (1 _ e)a 



r(e + v + \)T(e + v- A) 



-ea(x-L) + T(l + 2^)r(-l - 2e) e (l+ e )o(x-£) 



r(l + e + 1/ + A)r(l + e + i/ - A) T(-e + v + A)r(-e + i/ - A) 

Upon denning as S, T, U, V respectively the various combinations of gamma functions appearing 
the above expressions the wave functions are written as: 



(x) « A' 



g e -ea(x+L) , rp e -(l-e)a(x+L) 



(51a) 
(51b) 



Now we have to match the two solutions in x = requiring continuity of the wave-function 
0l(O) = ^_r(0) and of its first derivative 4>' L (0) = 4>' R (0). This gives the homogeneous system: 



A' 



Se -eaL +Te -(l-e)«L 



17' 



A' 



-e5 e" 



(1 - e)Te 



-(l-e)aL 



D' 



-eUe 



Ue +eaL + Ve -(l+e)aL 

aL + (1 + e)y e - (1+e)aL 









which admits a solution only if its determinant is zero. This provides a condition for extracting 
the energy eigenvalue: 



y ' TU 2e + l 



0. 



(52) 



When Eq. 52 is satisfied the relation between A' and D' is found to be: 

T It- 1 



D' 



e 2eaL A' = -e 

V 2e + 1 C/ 



2eaL j^l 



and j4' is the usual normalization constant. The condition in Eq. (52) is a transcendental equation 



which can be solved numerically. We provide numerical examples of the bound states. As we are 



studying bound states we seek numerical solutions of Eq. (52) in the interval — W < E < m. Since 



J-(E) = is complex the energy eigenvalues E are found by solving numerically, for real solutions, 
the two (independent) equations Ke[T(E)] = and rmfJ^E')] = 0. Fig. [4] shows graphically the 
details of the numerical computations both for the position dependent mass and the constant mass 
cases. In table I we give the numerical results for the spectra of the position dependent mass case 
and that of the constant mass case, with the same values of the parameters (mo = 1;L = 2; a = 10; 
W = 2) in order to make a meaningful comparison between the two cases. For the case of the 
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FIG. 4: The energy spectrum is derived by solving numerically with respect to the real variable E, the 
transcendental equations Re[J r (i?)] = 0, l\n[F(E)] = 0. In the figure we plot Ke[J 7 (E)} (solid line) and 
ImfJ^-E)] (dotted line) as functions of the energy E for the position dependent mass case (left panel) and 
for the constant mass case (right panel). The eigenvalues (full disks) are given by the points on the E'-axis 
where the two curves cross. The corresponding numerical values are given in table m 

constant mass we have used the results of |24j . We observe that the in the position dependent 
mass case the number of bound states decreases relative to the constant mass. In Figures [5] and 
[6] we provide some example of the (normalized) wave- functions and probability densities both for 
constant mass case and position dependent mass. Also, comparing Figure [5] with Figure [6] one can 
infer that in the position dependent mass case the probability density is almost flat in the region 
inside the potential well as opposed to the constant mass case where, for the highest excited states 
it oscillates strongly. We note, in the constant mass case (right panel of Figjl]), an eigenvalue 
corresponding to E = —m = — 1 which merges with the negative continuum. This situation has 
previously been considered in the literature [231 [291 EQ] and has been referred to as super-criticality. 
Such super-critical states are also called half-bound states and are characterized by the fact that 
one of the spinor components (the upper, u\, or the lower, 112) are not strictly normalizable. We 
show in table II] the eigenvalue E = — 1 only because it is a solution of Eq. 52 





E 1 


E 2 


E 3 


E 4 


E 5 


m = const. 


-1 


-0.759003 


-0.273555 


+0.271144 


+0.788942 


m(x) 


-0.633251 


-0.00806737 


+0.605869 







TABLE I: Numerical values of the energy eigenvalues (discrete spectrum) for the bound states. The model 
parameters are: mo = 1, W = 2, a = 10 and L = 2. In the upper row we have the spectrum of the constant 
mass case, while in lower row we have the position dependent mass case. 
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FIG. 5: The case of position dependent mass. Upper panel: plot of the normalized wave function i?e[iti] and 
of the probability density for E = —0.633251 (ground state); Lower panel: the same, but for E = +0.605869 
(highest excited state). The model parameters are: mo = 1, W = 2, a = 10 and L = 2. 



We do not address further the issue of super-criticality within the position dependent mass case 
as it goes beyond the scope of the present work. 

Finally in Fig [7] we give a further example in the case of the position dependent mass when 
the potential well is deeper, W = 3, with the other parameters as in Figure [5| We note that in 
this case the highest level is close to the continuum (E = mo = 1) and indeed the wave function 
converges less rapidly and the probability density as well. We have computed for example that in 
this case the probability of the particle to be outside the potential well (|x| > 2) is: P ou tside ~ 0.57 
which is even greater of the probability to be inside (— L < x < L) is: Pi ns ide ~ 0.43. The less 
rapid convergence of the wave function is due to the fact that the coefficient that controls such 
behavior (in this case as x — > oo,0#(x) ~ e~ uax ) is v = \J E 2 — rri^/a, and therefore very small 
giving a wave function that vanishes much slower than those corresponding to the lowest lying 
bound states. 
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FIG. 6: The case of constant mass. Upper panel: plot of the normalized wave function Re[u{\ and probability 
density for E = —0.759003 (lowest energy bound state); Lower panel: the same but for E = +0.7888942 
(highest excited state). The model parameters are the same as those of Figure [5} mo = 1, W = 2, a = 10 
and L = 2. 





FIG. 7: (Position dependent mass case). Plot of the normalized wave function Re[u2\ (left panel) and of 
the probability density (right panel) for E = 0.97248 (highest excited state). The model parameters are the 
same as those of Figure [5] except for W which now is W = 3. 
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V. DISCUSSION AND CONCLUSIONS 



We have solved the scattering problem for the one-dimensional Dirac equation with the WS 
potential in the position-dependent mass formalism. We have set some conditions on the equation 
in order to keep the structure of the hypergeometric equation which give a suitable mass function. 
These conditions provide a first order differential equation which can be solved exactly. With 
the physical requirement that the mass function at infinity goes to a constant mass mo, specifies 



completely the mass function. Once the mass function has been found, c.f. Eq. 25 , we have followed 
the same technique employed in ref. (23j solving the equation in the negative and positive region 
separately and giving the solution in terms of the hypergeometric function 2-^1- For the scattering 
problem ordinary boundary conditions at infinity are imposed and then the match at x = allows 
to specify all unknowns up to a normalization constant. 

We note that our method of solving the Dirac equation for a particular case of effective position 
dependent mass function has the drawback of providing a mass function which does not interpolate 
smoothly with the constant mass case. In other words our mass function does not contain a 
parameter such that when set to zero reduces the mass function to the constant case (mo). Further 
studies in this direction should be pursued in order to overcome such difficulties. 

We have obtained analytical expressions for the transmission and reflection coefficients, and 
we explicitly verified that unitarity (R + T = 1) is preserved in the PDM case. We have also 
studied the bound states, i.e. the discrete spectrum of the WS potential well with the effective 
position dependent mass, finding an exact analytical condition for the energy eigenvalues (in the 
form of a transcendental equation which needs to solved numerically). We have provided an 
explicit numerical example finding the eigenvalues and the wave function for a specific choice of 
the parameters. 

Our approach offers one of the few examples where the Dirac equation is solved exactly in the 
position dependent mass case and in an external potential. To the best of our knowledge the 
only other example is reported in ref. |18j where the three-dimensional Dirac equation is solved 
for the PDM case in the Coulomb external field. We note a similarity between the present work 
and that reported in |18| . In both cases the mass function for which an exact solution is found 
shares similarities with the external potential. In [18] the spherically symmetric mass function for 
which the problem is solved is m(r) = 1 + fi\ 2 jr where, in atomic units, mo = h = 1 and A is the 



Compton wavelength. Our mass function, c.f. Eq. 25 , is also certainly related to the shape of the 
Woods-Saxon external potential. 



22 



Acknowledgments 

This work has grown out of the diploma thesis of S.B. presented at the University of Perugia 
in September 2009. A. A. acknowledges warm hospitality from the Physics Department of the 
University of Perugia and INFN - Istituto Nazionale di Fisica Nucleare - Sezione di Perugia. 



[1] G. Bastard, Wave Mechanics Applied to Semiconductor Hetero- structures (EDP Sciences, Les Editions 
de Physique, Les Ulis, France, 1992). 

[2] O. von Roos, Phys. Rev. B 27, 7547 (1983). 

[3] T. Gora and F. Williams, Phys. Rev. 177, 1179 (1969). 

[4] O. von Roos and H. Mavromatis, Phys. Rev. B 31, 2294 (1985). 

[5] J. Thomsen, G. T. Eincvoll, and P. C. Hemmer, Phys. Rev. B 39, 12783 (1989). 

[6] J.-M. Levy-Leblond, Phys. Rev. A 52, 1845 (1995). 

[7] A. D. Alhaidari, Phys. Rev. A75, 042707 (2007), quant-ph/0703013. 

[8] A. de Souza Dutra and C.-S. Jia, Physics Letters A 352, 484 (2006), ISSN 0375-9601. 

[9] A. D. Alhaidari, H. Bahlouli, A. Al-Hasan, and M. S. Abdelmonem, Phys. Rev. A 75, 062711 (2007). 

[10] I. O. Vakarchuk, Journal of Physics A: Mathematical and General 38, 4727 (2005). 

[11] S. M. Ikhdair and R. Sever (2010), arXiv:1001.4327[math-ph]. 

[12] S. M. Ikhdair and R. Sever (2010), 1001.3943. 

[13] H. Egrifes and R. Sever, Int. J. Thcor. Phys. 46, 935 (2007). 

[14] C.-S. Jia, T. Chen, and L.-G. Cui, Physics Letters A 373, 1621 (2009), ISSN 0375-9601. 

[15] C.-S. Jia and A. de Souza Dutra, Annals of Physics 323, 566 (2008), ISSN 0003-4916. 

[16] L. Dekar, L. Chetouani, and T. F. Hammann, Journal of Mathematical Physics 39, 2551 (1998). 

[17] X.-L. Peng, J.-Y. Liu, and C.-S. Jia, Physics Letters A 352, 478 (2006), ISSN 0375-9601. 

[18] A. Alhaidari, Physics Letters A 322, 72 (2004). 

[19] K. S. Novoselov, A. K. Gchn, S. V. Morozov, D. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grigorieva, 
and A. A. Firsov, Science 306, 666 (2004), http://www.sciencemag.org/cgi/reprint/306/5696/666.pdf, 



URL http: //www. sciencemag. org/cgi/content/abstract/306/5696/666 
[20] K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, M. I. Katsnelson, I. V. Grigorieva, S. V. 



Dubonos, and A. A. Firsov, Nature 438, 197 (2005), URL doi : 10 . 1038/nature04233 



[21] A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, Rev. Mod. Phys. 81, 
109 (2009). 

[22] M. Li, S. Lee, and T. Kang, Current Applied Physics 9, 769 (2009), ISSN 1567-1739. 

[23] M. Setare and D. Jahani, Physica B: Condensed Matter 405, 1433 (2010), ISSN 0921-4526. 

[24] P. Kennedy, J. Phys. A35, 689 (2002), hep-th/0107170. 



23 



[25] M. I. Katsnelson, K. S. Novoselov, and A. K. Geim, NATURE PHYS. 2, 620 (2006), URL doi: 
|10.1038/nphys384) 

[26] M. Abramowitz and I. A. Stegun, Handbook of Matehematical Functions with Formulas, Graphs, and 

Mathematical Tables (Springer, 1992). 
[27] J. N. Ginocchio, Phys. Rev. Lett. 78, 436 (1997). 
[28] J. N. Ginocchio, Physics Reports 414, 165 (2005). 

[29] N. Dombey, P. Kennedy, and A. Calogeracos, Phys. Rev. Lett. 85, 1787 (2000). 
[30] A. Calogeracos and N. Dombey, Phys. Rev. Lett. 93, 180405 (2004). 



24 



